#include "header.h"

 static double udata[2004];
 static double vdata[2004];

//------------------------------------------------------------//

void initg(struct grid* in, int Np, double a, double dr){
  //Initializes one grid

  in->Np = Np;
  in->a  = a;
  in->dr = dr;
  in->b  = a+(Np)*dr; 
  in->rgrid   = (double*)malloc((Np+4)*sizeof(double)); 

  in->ugrid  = (double*)malloc((Np+4)*sizeof(double)); 
  in->vgrid  = (double*)malloc((Np+4)*sizeof(double)); 

  in->dgrid   = (double*)malloc((Np+4)*sizeof(double));  
	
  in->der1   = (double*)malloc((Np+4)*sizeof(double));  
  in->der1im   = (double*)malloc((Np+4)*sizeof(double));  

  in->der2   = (double*)malloc((Np+4)*sizeof(double)); 
  in->der2im   = (double*)malloc((Np+4)*sizeof(double)); 

  in->rhs     = (double*)malloc((Np+4)*sizeof(double)); 
  in->rhsim     = (double*)malloc((Np+4)*sizeof(double));

  in->tmp    = (double*)malloc((Np+4)*sizeof(double)); 
  in->tmpim    = (double*)malloc((Np+4)*sizeof(double));

  in->tmp2   = (double*)malloc((Np+4)*sizeof(double));
  in->tmp2im   = (double*)malloc((Np+4)*sizeof(double));

  printf("Refinement Success \n");

  in->rpointer=NULL;
  in->lpointer=NULL;

};

//------------------------------------------------------------//

void fillg(struct grid *in, double A, double sigma, int Gaussian){
  //fills one grid. A - amplitude of initial condition

  int  j;
  double a, dr;
  
  a  = in->a;  		// left endpoint
  dr = in->dr; 		// spatial step
  for (j=0; j<(in->Np)+4; j++){
    in->rgrid[j] = a + (j-2)*dr;  
    switch (Gaussian) {
      case 1:
	in->ugrid[j] = A*exp( -0.5*sigma* pow( (in->rgrid[j]), 2) ); 
	break;
		// Gaussian initial condition 
      case 2:
	in->ugrid[j] = A*exp( - pow( (in->rgrid[j])/sigma , 6) );	
	break;
		// Super-Gaussian IC in m
      case 3:
	if ( 1 - pow( (in->rgrid[j])/pow(25*A*A*sigma*sigma/3, 0.2) , 2) > 0 )
	{
	  in->ugrid[j] = ( A/pow(25*A*A*sigma*sigma/3, 0.3) )*sqrt(1 - pow( (in->rgrid[j])/pow(25*A*A*sigma*sigma/3, 0.2) , 2) )*exp(-0.2i*pow( in->rgrid[j], 2)/sigma );	
	} else {
	  in->ugrid[j] = 0;	
	}
	break;
		// Super-Gaussian IC in m

    }
  }
};

//------------------------------------------------------------//

void init_p(int Ngrids, struct grid* in){
  //initializes pointers

  int j;

  for (j=1; j<Ngrids-1; j++){
    in[j].rpointer = &in[j+1];
    in[j].lpointer = &in[j-1];
  }
  if (Ngrids!=1) {
    in[0].rpointer = &in[1];
    in[Ngrids-1].lpointer = &in[Ngrids-2];
  } 
}

//------------------------------------------------------------//

void glue(struct grid* in){
  //fills in ghost points on all grids

  struct grid *temp = in;
  struct grid *right, *left;

  while ((temp->lpointer)!=NULL) {temp = temp->lpointer;}

  while (temp != NULL){
 
    if ((temp->rpointer)!=NULL) {    
	right = temp->rpointer;
	temp->ugrid[temp->Np+2] = right->ugrid[2];
	temp->vgrid[temp->Np+2] = right->vgrid[2];

	temp->ugrid[temp->Np+3] = (5.*right->ugrid[2] + 15.*right->ugrid[3]
				 -5.*right->ugrid[4] + right->ugrid[5])/16.;
	temp->vgrid[temp->Np+3] = (5.*right->vgrid[2] + 15.*right->vgrid[3]
				 -5.*right->vgrid[4] + right->vgrid[5])/16.;


	temp->der1[temp->Np+2] = right->der1[2];
	temp->der1im[temp->Np+2] = right->der1im[2];

	temp->der1[temp->Np+3] = (5.*right->der1[2] + 15.*right->der1[3]
				 -5.*right->der1[4] + right->der1[5])/16.;
	temp->der1im[temp->Np+3] = (5.*right->der1im[2] + 15.*right->der1im[3]
				 -5.*right->der1im[4] + right->der1im[5])/16.;

	temp->der2[temp->Np+2] = right->der2[2];
	temp->der2im[temp->Np+2] = right->der2im[2];

	temp->der2[temp->Np+3] = (5.*right->der2[2] + 15.*right->der2[3]
				    -5.*right->der2[4] + right->der2[5])/16.;
	temp->der2im[temp->Np+3] = (5.*right->der2im[2] + 15.*right->der2im[3]
				    -5.*right->der2im[4] + right->der2im[5])/16.;


    } else {        
      temp->ugrid[temp->Np+2] = 2*temp->ugrid[temp->Np+1] - temp->ugrid[temp->Np];
      temp->vgrid[temp->Np+2] = 2*temp->vgrid[temp->Np+1] - temp->vgrid[temp->Np];

      temp->der1[temp->Np+2]  = 2*temp->der1[temp->Np+1] - temp->der1[temp->Np];
      temp->der1im[temp->Np+2]  = 2*temp->der1im[temp->Np+1] - temp->der1im[temp->Np];

      temp->der2[temp->Np+2]  = 2*temp->der2[temp->Np+1] - temp->der2[temp->Np];
      temp->der2im[temp->Np+2]  = 2*temp->der2im[temp->Np+1] - temp->der2im[temp->Np];
/*4.* self.mgrid[self.Np+1] - 6.*self.mgrid[self.Np]  
				+ 4.*self.mgrid[self.Np-1]  -  1.*self.mgrid[self.Np-2];*/

      temp->ugrid[temp->Np+3] = 2*temp->ugrid[temp->Np+1] - temp->ugrid[temp->Np-1];
      temp->vgrid[temp->Np+3] = 2*temp->vgrid[temp->Np+1] - temp->vgrid[temp->Np-1];

      temp->der1[temp->Np+3]  = 2*temp->der1[temp->Np+1] - temp->der1[temp->Np-1];
      temp->der1im[temp->Np+3]  = 2*temp->der1im[temp->Np+1] - temp->der1im[temp->Np-1];

      temp->der2[temp->Np+3]  = 2*temp->der2[temp->Np+1] - temp->der2[temp->Np-1];
      temp->der2im[temp->Np+3]  = 2*temp->der2im[temp->Np+1] - temp->der2im[temp->Np-1];

/* 8.*self.mgrid[self.Np+2] - 22.*self.mgrid[self.Np+1]
				 + 28.*self.mgrid[self.Np] - 17.*self.mgrid[self.Np-1]  + 4*self.mgrid[self.Np-2];*/
    }
    if ((temp->lpointer)!=NULL) {
	left  = temp->lpointer;
	temp->ugrid[0] = left->ugrid[left->Np-2];
	temp->vgrid[0] = left->vgrid[left->Np-2];

	temp->der1[0] = left->der1[left->Np-2];
	temp->der1im[0] = left->der1im[left->Np-2];

	temp->der2[0] = left->der2[left->Np-2];
	temp->der2im[0] = left->der2im[left->Np-2];

	temp->ugrid[1] = left->ugrid[left->Np]; 
	temp->vgrid[1] = left->vgrid[left->Np]; 

	temp->der1[1] = left->der1[left->Np];
	temp->der1im[1] = left->der1im[left->Np];

	temp->der2[1] = left->der2[left->Np];
	temp->der2im[1] = left->der2im[left->Np];

      } else {  
	  temp->ugrid[0] = temp->ugrid[4];
	  temp->vgrid[0] = temp->vgrid[4];

	  temp->der1[0] = -temp->der1[4]; 
	  temp->der1im[0] = -temp->der1im[4]; 

	  temp->der2[0] =  temp->der2[4];
	  temp->der2im[0] =  temp->der2im[4];

	  temp->ugrid[1] = temp->ugrid[3];
	  temp->vgrid[1] = temp->vgrid[3];

	  temp->der1[1] = -temp->der1[3]; 
	  temp->der1im[1] = -temp->der1im[3]; 

	  temp->der2[1] =  temp->der2[3];
	  temp->der2im[1] =  temp->der2im[3];
	
      }
    temp = (temp->rpointer);
   }
}

//------------------------------------------------------------//

void d1r(struct grid* in){
  // computes 4-point stencil derivative in interior               
  int k;
  struct grid *temp = in;
   
  while ((temp->lpointer)!=NULL) {temp = temp->lpointer;}

  while (temp != NULL){

    for (k=2; k<((temp->Np)+2); k++){

      temp->der1[k] = (-temp->ugrid[k+2] + 8.*temp->ugrid[k+1]
			 - 8.*temp->ugrid[k-1] + temp->ugrid[k-2]);
      temp->der1im[k] = (-temp->vgrid[k+2] + 8.*temp->vgrid[k+1]
			 - 8.*temp->vgrid[k-1] + temp->vgrid[k-2]);

      temp->der1[k] = temp->der1[k]/(12.*temp->dr);

      temp->der1im[k] = temp->der1im[k]/(12.*temp->dr);
    }


    temp = (temp->rpointer);
  }
  glue(in);
}

//------------------------------------------------------------//

void d2r(struct grid* in){
  //computes 5-point stencil 2nd derivative in interior             
  int k;
  struct grid *temp = in;
   
  while ((temp->lpointer)!=NULL) {temp = temp->lpointer;}

  while (temp != NULL){

    for (k=2; k<((temp->Np)+2); k++){

      temp->der2[k] = (- temp->ugrid[k+2] + 16.*temp->ugrid[k+1] 
		       - 30.*temp->ugrid[k] + 16.*temp->ugrid[k-1] - temp->ugrid[k-2]);

      temp->der2im[k] = (- temp->vgrid[k+2] + 16.*temp->vgrid[k+1] 
		       - 30.*temp->vgrid[k] + 16.*temp->vgrid[k-1] - temp->vgrid[k-2]);

      temp->der2[k] = temp->der2[k]/(12.*temp->dr*temp->dr);

      temp->der2im[k] = temp->der2im[k]/(12.*temp->dr*temp->dr);
    }


    temp = (temp->rpointer); 
  }
  glue(in);
}


//------------------------------------------------------------//

//------------------------------------------------------------//
void rhs(double a, struct grid* in){
  // computes right hand-side on all grids                           
  struct grid *temp = in;
  struct grid *right;
  int k;
   
  //glue(in);
  while ((temp->lpointer)!=NULL) {temp = temp->lpointer;}
  right = temp;
  d2r(temp);  
  d1r(temp);  

  temp = right;
  while (temp != NULL){

     for (k=2; k<((temp->Np)+2); k++){

       temp->rhs[k] =  - temp->der2im[k] - 2.*((temp->der1im[k])/(temp->rgrid[k])) 
	 - ( pow(temp->ugrid[k], 2) + pow(temp->vgrid[k], 2) )*(temp->vgrid[k]) ;

       temp->rhsim[k] =  temp->der2[k] + 2.*((temp->der1[k])/(temp->rgrid[k])) 
	 + ( pow(temp->ugrid[k], 2) + pow(temp->vgrid[k], 2) )*(temp->ugrid[k]) ;
     } 
     temp = (temp->rpointer);
  }
  temp = in;
  while ((temp->rpointer)!=NULL) {temp = temp->rpointer;}
  temp->rhs[temp->Np+2] = 0;
  temp->rhsim[temp->Np+2] = 0;

  temp->rhs[temp->Np+3] = 0;
  temp->rhsim[temp->Np+3] = 0;

  while ((temp->lpointer)!=NULL) {temp = temp->lpointer;}
  temp->rhs[2] =  - 3.*(temp->der2im[2]) - ( pow(temp->ugrid[2], 2) + pow(temp->vgrid[2], 2) )*temp->vgrid[2] ;
  temp->rhsim[2] =  3.*(temp->der2[2]) + ( pow(temp->ugrid[2], 2) + pow(temp->vgrid[2], 2) )*temp->ugrid[2];
}

//------------------------------------------------------------//

double Nparticle(struct grid* in){
  struct grid *temp = in;
  double Npart = 0;
  int i = 0;

  while ((temp->lpointer)!=NULL) {temp = temp->lpointer;}

  while (temp != NULL ) {
    for ( i = 2; i < (temp->Np) + 3; i++ ) {
       Npart = Npart + 0.5*(temp->dr)*(temp->rgrid[i])*(temp->rgrid[i])*(pow(temp->ugrid[i],2) + pow(temp->vgrid[i],2));
       Npart = Npart + 0.5*(temp->dr)*(temp->rgrid[i])*(temp->rgrid[i+1])*(pow(temp->ugrid[i+1],2) + pow(temp->vgrid[i+1],2));
	
    }
    temp =  temp->rpointer;	
  } 
 return Npart;	
	
}

//------------------------------------------------------------//
/*void density(struct grid* in){
  // fills in array of densities                                      
  struct grid *temp = in;
  struct grid *left;
  int i;
   
  while ((temp->lpointer)!=NULL) {temp = temp->lpointer;}
  
  left = temp;

  d1r(temp);

  temp = left;

  while (temp != NULL){

    for (i=2; i<(temp->Np+2); i++){
      temp->dgrid[i] = temp->der1[i]/(temp->rgrid[i]);
    }
    
    temp = (temp->rpointer);
  } 

  left->dgrid[2] = (-left->mgrid[4] + 16.*left->mgrid[3] 
		    -30.*left->mgrid[2] + 16.*left->mgrid[1] - left->mgrid[0]);

  left->dgrid[2] = left->dgrid[2]/(12.*left->dr*left->dr);

  }*/

//------------------------------------------------------------//

void clearall(struct grid *in){
  // calls free for all grids
  struct grid *temp = in;

  while (temp != NULL){
    free(temp->rgrid);
    free(temp->ugrid); 
    free(temp->vgrid);
    free(temp->dgrid);
    free(temp->der1); 
    free(temp->der1im);
    free(temp->der2); 
    free(temp->der2im);
    free(temp->rhs); 
    free(temp->rhsim);

    temp = (temp->rpointer);
  }
}

//----------------------------------------------------
void addgrid(struct grid* in, int *Ngrids){
  // links one grid to array of grids, does interpolation
  int i;
  struct grid* temp = in;
  struct grid* right; 

  while ((temp->lpointer)!=NULL ){
    temp = (temp->lpointer);
  }

  initg(&in[*Ngrids], temp->Np, 0., 0.5*temp->dr );
  fillg(&in[*Ngrids], 0., 0, 1.);

  in[*Ngrids].rpointer = temp;

  for (i = 1; i < 2 + (temp->Np)/2; i++ ){

    in[*Ngrids].ugrid[2*i] = temp->ugrid[i+1];
    in[*Ngrids].vgrid[2*i] = temp->vgrid[i+1];

    in[*Ngrids].ugrid[2*i + 1] = (3*temp->ugrid[i-1] - 25*temp->ugrid[i] + 150*temp->ugrid[i+1]
				  + 150*temp->ugrid[i+2] - 25*temp->ugrid[i+3] + 3*temp->ugrid[i+4])/256.;

    in[*Ngrids].vgrid[2*i + 1] = (3*temp->vgrid[i-1] - 25*temp->vgrid[i] + 150*temp->vgrid[i+1]
				  + 150*temp->vgrid[i+2] - 25*temp->vgrid[i+3] + 3*temp->vgrid[i+4])/256.;
  }

  
  in[*Ngrids].ugrid[1] = in[*Ngrids].ugrid[3];
  in[*Ngrids].ugrid[0] = in[*Ngrids].ugrid[4];

  in[*Ngrids].vgrid[1] = in[*Ngrids].vgrid[3];
  in[*Ngrids].vgrid[0] = in[*Ngrids].vgrid[4];

  for (i = 0; i < 4 + (temp->Np)/2; i++ ){

    udata[i] = temp->ugrid[i + (temp->Np)/2];
    vdata[i] = temp->vgrid[i + (temp->Np)/2];
  }
  right = temp->rpointer;
  
  temp->Np = (temp->Np)/2; 
  temp->a  = in[*Ngrids].b;

  initg(temp, temp->Np, temp->a, temp->dr);
  temp->lpointer = &in[*Ngrids];
  temp->rpointer = right;
  fillg(temp, 0., 0, 0.);
  
  //
  for (i = 0; i < 4 + temp->Np; i++ ){
    temp->ugrid[i] = udata[i];
    temp->vgrid[i] = vdata[i];
  }

  glue(in);
  ++*Ngrids;
}

//-------------------------------------------------------------------
void RK4(struct grid* in, double* dt, double a){
  struct grid* temp = in;
  struct grid* left;
  int i;

  while ((temp->lpointer)!=NULL) {temp = temp->lpointer;}
  left = temp;
  rhs(a, temp);

  while (temp!=NULL){
    for (i = 0 ; i < temp->Np+3 ; i++){
      temp->tmp[i]   = temp->ugrid[i] + (*dt/(6.)) * temp->rhs[i];
      temp->tmpim[i]   = temp->vgrid[i] + (*dt/(6.)) * temp->rhsim[i];

      temp->tmp2[i]  = temp->ugrid[i];
      temp->tmp2im[i]  = temp->vgrid[i];

      temp->ugrid[i] = temp->ugrid[i] + (*dt/(2.)) * temp->rhs[i];
      temp->vgrid[i] = temp->vgrid[i] + (*dt/(2.)) * temp->rhsim[i];
    }
    temp = temp->rpointer;
  }  
  temp = left;
  glue(temp);
  temp = left;
  rhs(a, temp);
  temp = left;

  while (temp!=NULL){
    for (i = 0 ; i < temp->Np+3 ; i++){
      temp->tmp[i]   = temp->ugrid[i]  + (*dt/(3.)) * temp->rhs[i];
      temp->tmpim[i]   = temp->vgrid[i]  + (*dt/(3.)) * temp->rhsim[i];

      temp->ugrid[i] = temp->tmp2[i]   + (*dt/(2.)) * temp->rhs[i];
      temp->vgrid[i] = temp->tmp2im[i]   + (*dt/(2.)) * temp->rhsim[i];
    }
    temp = temp->rpointer;
  }  
  temp = left;
  glue(temp);  
  temp = left;
  rhs(a, temp);
  temp = left;

  while (temp!=NULL){
    for (i = 0 ; i < temp->Np+3 ; i++){
      temp->tmp[i]   = temp->ugrid[i]  + (*dt/(3.)) * temp->rhs[i];
      temp->tmpim[i]   = temp->vgrid[i]  + (*dt/(3.)) * temp->rhsim[i];

      temp->ugrid[i] = temp->tmp2[i]   + (*dt) * temp->rhs[i];
      temp->vgrid[i] = temp->tmp2im[i]   + (*dt) * temp->rhsim[i];
    }
    temp = temp->rpointer;
  }  
  temp = left;
  glue(temp);
  temp = left;
  rhs(a, temp);
  temp = left;
  while (temp!=NULL){
    for (i = 0 ; i < temp->Np+3 ; i++){
      temp->ugrid[i] = temp->tmp[i]  +  (*dt/(6.)) * temp->rhs[i];
      temp->vgrid[i] = temp->tmpim[i]  +  (*dt/(6.)) * temp->rhsim[i];
    }
    temp = temp->rpointer;
  }  
  temp = left;
  glue(temp);
  temp = left;

}

//--------------------------------------------------//


void startup(struct grid *in, struct INPdata *input, int *ac, char ***av){
  // setting up run
  char filename[80];
  FILE* fh;
  double a, dr;
  int i;

  if (*ac != 2) {
   printf( "\n\t Usage: %s run_name\n\n", *av[0]);
   printf( "\n\t\t %s run_name binary_file\n\n", *av[0]);
   exit( 1 );
  }
  printf("\nhere\n  " );

  strcpy( filename, input->runname );
  printf("\nhere\n  " );
  strcat( filename, ".inp" );
  printf("%s" ,  input->runname);

  fh = fopen( filename , "rt" );
  if ( fh == NULL ) {
    printf ( "\n  File \"%s\" not found.\n\n", filename);
    exit( 1 );
  }
  read_p( fh, input );
  fclose( fh );

  a = 0.;                                                       //  grid  0
  dr = input->dr;                                              // spacing of grid 0

  for (i = 0; i <  (input->Ngrids); i++){

    initg( &in[i] , (input->Npoints[i]) , a, dr );                   // setup grid spacing
    init_p( input->Ngrids , in );                                  // setup pointers
    fillg( &in[i], input->A, input->sigma, input->Gaussian );      // fill initial data 

    a  = in[i].b;                                                  //  next grid 
    dr = 2.*dr;                                                  //  next grid
  }

  if ( input->Gaussian == 3 ) {
    sgaussIC( in, &(input->Ngrids), &(input->A), &(input->sigma));
    glue( in );
  }
  //density( in );
  
}

